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Abstract 

Universal properties of the zero temperature superconductor-insulator tran- 
sition in two-dimensional amorphous films are studied by extensive Monte 
Carlo simulations of bosons in a disordered medium. We report results for 
both short-range and long-range Coulomb interactions for several different 
points in parameter space. In all cases we observe a transition from a super- 
conducting phase to an insulating Bose glass phase. Prom finite-size scaling 
of our Monte Carlo data we determine the universal conductivity a* and the 
critical exponents at the transition. The result a* = (0.55 ± 0.06)(2e) 2 //i for 
bosons with long-range Coulomb interaction is roughly consistent with exper- 
iments reported so far. We also find a* = (0.14±0.03)(2e) 2 //i for bosons with 
short-range interactions. 

PACS numbers: 74.65.+n, 74.70.Mq, 74.75.+t 
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I. INTRODUCTION 



From the work of Abrahams, Anderson, Licciardello, and RamakrishnanB it is known 
that no true metallic behavior can be observed for non-interacting electrons at T = 
in two dimensions, since all states will be localized by arbitrarily weak disorder. When 
repulsive interactions are turned on the situation is less clear but the general beliei is that 
a metallic phase still should be absent at T = in the presence of disorder, although we 
know of no rigorous proof of this. However, in the presence of attractive interactions, a 
superconducting phase is expected!, both at T = and finite T, even for a finite amount 
of disorder, because disorder is irrelevant^ at the finite temperature transition, which is of 
the Kosterlitz-Thoulessi type discussed below. The onset of superconductivity at T = 
is presumed then, in d = 2, to be directly from the insulating phase with no intervening 
metallic phase. One should therefore in principle be able to observe a direct insulator- 
superconductor transition at zero temperature in two dimensions as a function of disorder 
and/or interaction strength. The main topic of this paper is to analyze such a transition 
and extract its universal features. 

Dimensionality and divergent length scales play an important role in continuous phase 
transitions. The diverging correlation length scale implies that many microscopic details are 
irrelevant. Furthermore physical quantities containing dimensions of length to some non-zero 
power typically diverge or vanish at the critical point. Two dimensions is special in that the 
conductivity contains no length scale units, i.e. the conductance per square is the same as 
the conductivity. Hence, right at the T = quantum critical point, the conductivity is not 
only finite and nonzero but also universal^, even though it is zero in the insulating phase 
and infinite in the superconducting phase. This view differs from that of previous workll 
which parameterized the transition in terms of the normal state resistivity. The calculation 
of this universal conductivity is one of the main goals of the present paper. A short account 
on some of our results has already been publishedi. 

A schematic phase diagram is shown in Fig. [I] as a function of temperature, T, and 
disorder, A. At zero temperature, a critical amount of disorder, A c , separates the supercon- 
ducting from the insulating phase. Even at finite temperatures the superconducting phase 
persists whereas a truly insulating phase only exists at T = 0, because, at finite T, electrons 
can be inelastically scattered from one localized state to another, and hence conduct. This 
insulating phase, consisting of localized electron pairs, can then be described, close to the 
critical point, as a Bose-condensed fluid of vortices. The universality class of the transition 
should therefore be that of the superconductor to Bose glassS. 

Let us first discuss the nature of the transition at finite temperatures indicated by the 
solid line in Fig. [TJ. The finite temperature transition should have many similarities with the 
2D XY transition at which logarithmically interacting vortices unbindi. However, Pearl0 
showed that in a superconducting film vortex pairs only have logarithmic interactions out 
to a distance = 2X 2 /d beyond which the interaction energy falls off as 1/r. Here A is 
the bulk penetration depth, d the film thickness, and the screening length for magnetic 
fields. Due to this cutoff, the energy required to create a vortex is always finite and no sharp 
transition should exist. However, according to the Kosterlitz-Thouless theory!, the value 
of A ± (T C ) is given exactly bytH3 A ± (T C ) = (f) 2 /(lQir 2 k B T c ), where O = hc/2e is the flux 
quantum. Numerically A±(T C ) = 2/T c where T c is in Kelvin and A±(T C ) is in centimeters. 
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Thus, Aj_ is so large at T c that rounding of the transition due to the presence of free vortices 
below T c is almost certainly unobservable. In fact, rounding due to finite-size effects is 
probably more important. 

The vortex unbinding transition at T c is driven by fluctuations in the phase of the super- 
conducting order parameter. At a higher temperature, T c o, fluctuations in the amplitude of 
the order parameter will become important and a crossover to a regime dominated by para- 
conduct ivityEH will occur. T c0 is indicated by the dashed line in Fig. [IJ. For T c < T <C T c0 the 
presence of free vortices destroys the characteristic global properties of the superconducting 
phase. Nevertheless, a local order parameter still exists between T c and T c0 . The presence 
of free vortices leads to a finite conductivity of the formlll a v ~ 0.37cr n (£/£ c ) 2 , where a n is 
the conductivity of the normal state electrons, £ c the core size of a vortex and £, a typical 
distance between free vortices, is the Kosterlitz-Thoulessll correlation length which diverges 
exponentially at T c . The exponential tail in the resistivity caused by the presence of free 
vortices between T c and T c0 has been observed experimentally0'0 in superconducting films 
with high normal state resistivity. In these experiments the "mean field onset temperature" , 
T c q is determined by fitting the resistance to an Aslamazov-LarkinllS form and T c o — T c is 
found to be of the order of half a Kelvin. In the dirty limit, Beasley et al.B-3 derive the rela- 
tionship t c = (T c0 — T c )/T c ~ 0.17e 2 /cr n h so, for films with a relatively high sheet resistance, 
t c can be appreciable. For a review of the finite temperature transition we refer the reader 
to Refs. [HJ More recently some evidence for a vortex-antivortex unbinding transition in 
superconducting niobium filmsS with R n = has been found. However, the difference 

between T c and T c q is very small in the clean limit which makes the Kosterlitz-Thouless 
behavior difficult to observe. 

The superconducting order parameter is a complex scalar, described by both a magnitude 
and a phase. Our key basic assumption is that universal properties at superconductor- 
insulator transition are determined only by phase fluctuations, as outlined above, and that 
the magnitude of the order parameter, and therefore of the gap in the fermionic energy 
spectrum, remains finite at the critical point. We thus assume that when disorder drives 
T c to zero, T c0 remains nonzero. If the transition is approached from the insulating side, a 
local order parameter appears before the onset of global phase coherence at A c . Implicit in 
our assumption is that Cooper pairs, and thus a gap to fermionic excitations, persists into 
the insulating phase, even though superconductivity is destroyed by phase fluctuations. On 
the scale of a diverging phas e-correlation length, £, the individual Cooper pairs will look like 
point particles. The fermionic degrees of freedom should therefore be highly suppressed at 
the critical point and an approximate description in terms of point-like bosons should be 
valid. It is possible that strong disorder destroys the local fermionic gap at a finite density 
of points, but, provided that the Fermi degrees of freedom are localized, they may still be 
dynamically irrelevant and our model applicable. 

Since we shall be concerned with the transition at A = A c , T = 0, vortices in the 
system will not be excited thermally, but there will be vortices present created by quantum 
fluctuations^. We therefore need to treat the vortices as quantum mechanical objects and 
one might expect the transition at A c to be described by a (2+l)D XY modelil, where the 
extra dimension arises because we are considering a T=0 quantum phase transition, see e.g. 
Ref. [K]. However as we shall see in Section [H], although the physics is indeed described 
by a (2+l)D system, its symmetry is, in general, not that of the XY model. It is also 
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worth noting that the vortex mobility, (2e 2 /7r/i 2 )£ 2 i? n lll, is significantly augmented in dirty 
superconducting films. Hence, at A c the vortices should be seen as fairly light objects that 
move rather freely. At still higher disorder the Bose glass phase should cross over into an 
Fermi glass when the individual electrons constituting the bosons become localized. This 
behavior may have already been observed in a magnetic fieloS. 

Recent experiments seem to confirm that a direct insulator-superconductor transition 
indeed does take place at zero temperature in many materials. Haviland et al.@ and Liu et 
al.il have performed experiments on Bi films grown in situ. The experimental technique is 



described in Ref. |26|. These films are believed to be truly amorphous on an atomic scale. 
The authors report a critical d.c. resistivity, R^, very close to Rq, where 

R Q = h/Ae 2 w 6453fi . (1.1) 

Furthermore, experiments performed on DyBaCuO films§§ and NdCeCuOH show clear 
evidence for a direct superconductor-insulator transition. The reported critical resistivity 
seems in this case to be somewhat higher, around 10 kO or 1.5 Rq EE. Lee and Ket- 
tersonil have presented results from experiments on MoC films again showing very clear 
evidence for a superconductor-insulator transition occurring at zero temperature, but with 
i?n slightly lower, in the range 2.8-3.5 kf2 ~ 0.5-Rq. Furthermore, experiments performed 
on Josephson junction arrays@'0, which are believed to be in the same universality class, 
also seem to support the picture of a superconductor-insulator transition. The existence 
of a superconductor-insulator transition in two-dimensional films at zero temperature thus 
seems well established but evidence for the universality of the critical resistivity remains 
weak. It is not clear, however, whether all the experiments are in the critical region. In 
order to establish that a given experiment is actually probing the critical regime, one must 
show scaling of the resistivity data. This has been done successfully by Hebard and Paala- 
ner£§o for the field-tuned transition and partially successfully by the Minnesota grou pi. 
However it is likely that most measurements to date have failed to probe the critical regime, 
and further experiments at even lower temperatures are expected to give better agreement 
among the different estimates of R^. 

The situation concerning the relevance of a bosonic picture seems less clear. Hebard and 
PaalaneniE have reported results on amorphous InO films in a magnetic field, supporting 
the existence of Cooper pairs in the insulating phase. For the B = transition, Hebard and 
PaalanenS have presented clear transport evidence that T c q remains finite as T c is driven 
to zero. On the other hand, direct tunneling measurements by Dynes et al.0 and Valles et 
al.il on homogeneously disordered Pb films shows that the gap goes to zero at the critical 
point. There seems, however, to be a general agreement that a local superconducting order 
parameter exists prior to the transition in granular films and in Josephson junction arrays 
where the individual grains become superconducting above T c . It is possible that tunneling 
experiments tend to emphasize regions of the samples containing quasilocalized fermion 
states below the gap which are necessary to achieve tunneling. 

A number of theoretical and numerical studies of the superconductor-insulator transition 
have been performed. GoldH studied the impurity induced insulating transition in the 
interacting Bose gas. Giamarchi and SchulzS considered the one-dimensional electron gas 
with attractive interactions in the presence of disorder. They found a transition to a localized 
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phase in the same universality class as that of repulsively interacting bosons in a random 
potential. This lends strong support to our assumption of the dirty boson universality class 
in the 2D case. Fisher et al.clTcall2l considered the boson Hubbard model and, through a 
scaling analysis, derived equations for the exponents governing the superconductor-insulator 
transition as well as the phase diagram for dirty bosons. A renormalization group approach 
was taken by Weichman et al@ who performed a double epsilon expansion for the dirty boson 
problem. Following the initial suggestion of a Bose glassphase in the disordered system and 
a Mott insulator in the clean system, Batrouni et al.@ and Krauth et al.@, showed, by 
quantum Monte Carlo simulations, the existence of Mott insulating phases in an interacting 
boson system without disorder, characterized by the exponents predicted by Fisher et al.H2l 
Subsequently, these authors considered the disordered case and evidence for a Bose glass 
was found.oa A Bose glass phase was also observed in a real space renormalization group 
study by Singh et al.0 The universal conductivity was first calculated by a 1/N expansion 
and Monte Carlo methods for the (2+l)D XY model by Cha et al.0 and Girvin et al.0H 
The universal conductivity for disordered bosons was then calculated by Rungell by exact 
diagonalization techniques on small lattices. Universal properties for a boson system in 
the presence of disorder both with and without long-range interactions were calculated by 
S0rensen et alS0 by Monte Carlo simulations, using a path integral representation which, 
effectively, only includes phase fluctuations in the Bose field. A universal conductivity was 
also recently found by Kampf et al.il in the boson Hubbard model including both phase and 
amplitude fluctuations. Two recent works have recently been published after the present 
work was finished. Batrouni et al.@ have calculated the universal conductivity by quantum 
Monte Carlo simulations directly on the boson Hubbard model, and Makivic et al.E3 have 
calculated the exponents and the universal conductivity using a hard-core boson model. The 
results of these last two papers differ from ours, and we shall comment on this in section 



VIII. 



Here we shall consider two forms of interaction between the bosons: short-range repulsive 
interaction and long-range Coulomb interaction. The model with short-range interactions is 
relevant to experiments on the onset of superfluidity in 4 He filmsS However, our present re- 
sults for the zero temperature transition in 2D are not directly applicable to He experiments 
in porous media such as Vycor or xerogelSH, since these experiments are mainly concerned 
with the 3D transition at finite temperatures. As stated above, the model with Coulomb 
interactions is expected to be in the correct universality class to describe the superconductor- 
insulator transition. However, the model and many of the results presented in this paper are 
applicable to other systems too. The world lines of the dirty boson model describe a gas of 
stringlike objects in a random medium. In addition to the superconductor-insulator transi- 
tion this model may also apply to other problems such as vortex lines in high-temperature 
superconductors with correlated pinning centersi^E, and polymer solutions. Our results for 
universal quantities might also be relevant for these problems. 

The organization of the paper is as follows. In Section [IJ we shall construct a form of the 
boson Hubbard model, including disorder and interactions, which is suitable for Monte Carlo 
simulations. Here we shall assume, as discussed above, that only bosonic degrees of freedom, 
i.e. complex order parameter fluctuations, are relevant at the superconductor-insulator tran- 
sition. In addition, to further simplify the numerical work, we effectively include only phase 
fluctuations of the bosons, amplitude fluctuations being neglected. Section |TJ describes the 
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scaling theory of the quantities that we are interested in. We discuss, in Section [TV], how we 
determine these quantities in the simulation, and we also treat the finite-size scaling tech- 
niques needed to extrapolate our results to infinite size. Section [V] describes our Monte Carlo 
methods, while Section [V| presents our results for short-range interactions and disorder, rel- 
evant to experiments on helium films. In Section |V11| long-range Coulomb interactions are 
included along with disorder. We believe that this model contains all ingredients necessary 
to make it relevant to experiments on the superconductor-insulator transition; i.e., that it 
is in the correct universality class. Our results are discussed in Section [VIII. 



II. THE MODEL 

In this section we introduce our basic model and via a sequence of transformations arrive 
at a form suitable for Monte Carlo simulation. As argued above it should be possible to 
describe the universal features of the superconductor-insulator transition in terms of boson 
physics. In this section we shall argue that the relevant starting point is the boson Hubbard 
model with a random local chemical potential (site energy). If only phase fluctuations are 
relevant we can map this model onto a dual Villain type model. We shall see that only in 
the absence of disorder and when there is an integer number of bosons per site, does this 
model belong to the same universality class as the (2+l)D XY model. 

In order to model the superconductor-insulator transition in terms of bosons, we must 
include an on-site repulsive interaction, otherwise all bosons would collapse into the lowest 
lying, highly localized state. The on-site repulsion term is the simplest possible way to model 
Coulomb repulsion. The correct treatment of the long-range part of the interaction will be 
discussed below. For simplicity we shall take an underlying square lattice of spatial size 
N = L x L. Changes in the symmetry of the lattice are not expected to modify the critical 
behavior of the model. We can then write down the boson Hubbard model in presence of 
disorderSHH: 



HbH = H + Hi (2.1) 



where 



* r r 

E x = -t (&A> + MJO • (2.2) 

<r,r'> 

Here U is the on-site repulsion, \i is the chemical potential, z the number of nearest neighbors, 
and v r represents the random on-site potential varying uniformly in space between —A and 
A. As usual, n r = &l& r is the number operator on site r. The hopping strength is given by 
t, and (r, r') denotes summation over pairs of nearest neighbors, each pair counted once. 

In the absence of disorder there is no insulating phase unless we fix the boson density at 
an integer value, no- Let us consider this case first. If we set $ r = |<i> r |e ?,e ' r and integrate out 
amplitude fluctuations the boson Hubbard model, Eq. (|2.2|) , becomes a model of coupled 
Josephson junctions,li3ELl 
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#jj = j E K - E * cos (^ - e r ,) , (2.3) 

(r,r'> 

where, in this representation, n r , which denotes the deviation of the boson number from no, 
runs from — oo to oo and so Eq. ( [2.3|) can only be quantitatively compared with the Hubbard 
model, Eq. (|2.2|), when no is very large, but is expected to be in the same universality class 



for arbitrary integer no- Note that t in Eq. (2J3) is 2no times the parameter t in Eq. ( [2.2| ). 
The phase operator, 8 r is canonically conjugate to n r so this version of the boson Hubbard 
model can be written in the angle representation as the quantum rotor modeSffl! 

(2-4) 

Let us write the partition function corresponding to H qr as 

Z = Tr exp[-/3(T + V)} , (2.5) 
where the kinetic energy of the rotors is 

U d 2 

(which corresponds to the potential energy of the bosons) and the potential energy of the 
rotors is 

V = - E tcos(6 r -6 T ,) . (2.7) 

<r,r'> 

We evaluate the trace in the partition function by writing a path integral over M time 
slices Tj between r = and r = (3: 



Z = Ti {exp[-/3(T + V)}/M} M 



= lim Tr{exp[-Ar T}exp[-ArV}} M , (2.8) 

M— >oo 

where hr is imaginary time and 

At = p/M (2.9) 

is the width of one time slice. Note that the limit At — > must be taken to correctly rep- 
resent the underlying quantum mechanics problem. Eq. (|2.8|) can be rewritten by inserting 
complete sets of states 



. Af-l 

IveU ({Q{Tj+i)}\ exp[— ArT] exp[— AtV]|{0(tj)}) , (2.10) 

3=0 



where \{0{Tj)}) is a coherent state in which site r has phase 9 r (T 3 ) at time Tj and the trace 
is enforced by periodic boundary conditions 
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Wtm)} = {6(t )} . 
The coherent states are eigenstates of the potential so 



(2.11) 



exp[-ArV] {0(7})}} = exp { Art £ cos frfa) - 6 r {r 3 )] ) {0(t})} 

<r,r') 



(2.12) 



where the sum is over all nearest neighbor spatial pairs, and hence Eq. fl2.10|) becomes 



M-l 



where 



and 



Z « fve n exp K x cos^Arj) ~ U^)) Tj , 

3=0 [ <r,r') I 



2i = ({^(r J+ i)}|e- A ^|{^(r,)}) 



(2.13) 



(2.14) 



K x = tAr . (2.15) 
Since the kinetic energies on different sites commute, we can consider each site separately: 

~AtU d 2 



exp 



2 del 



(2.16) 



Let Jj!" (rj) be the integer- valued angular momentum at r at time Tj. The corresponding 
state has wave function 



(O r (T j )\JZ(T j )) = e iJ Zto»'M 



(2.17) 



which is an eigenfunction of the kinetic energy. Inserting this complete set of states, we have 

ll^,(7- y+ i)|-A-(',»;^i>i- A/( 



Tj = EIl(^ + i)l^))exp {-^ [j;(r,)] 2 } (j;(r,)|^(r,)) , (2.18) 



and thus 



M-l 



[V9J2^P\K X J2 H cos [0^7}) - 0,(7})] 

{J} [ {r,r'> J=0 



x exp 



Art/ 



M-l 



E E tffo)]' 



j=0 



x exp { z£ E Prfa) - r (r i+1 )] 

r j=0 



(2.19) 



We can now proceed in two possible ways. We can either integrate out the angular variables 
{0} to obtain a statistical mechanics problem in the integer variables {J}, or we can sum 
over the {J} to obtain a classical (2+l)-dimensional XY model. Let us start with the latter. 
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Because Ar is small, the sum over the {</} is slowly convergent. We may remedy this 
by using the Poisson summation formula 

oo r+oo 

F{0) =Y,e-^ TUj2/2 e tje = >T / dJe 2mJm e- ATUj2 / 2 e lJd 



m=— 00 
00 



E \l^'^ {9 ' 2nm)2 ■ ( 2 - 2 0) 



m=— 00 



This periodic sequence of narrow Gaussians is the Villain approximation^ to the periodic 
function 

F(0) w e ^ cos(e) , (2.21) 
where we have dropped an irrelevant constant prefactor, and 

K T = — — . (2.22) 



Using this result in Eq. (|2.19 ) we finally arrive at the partition function of the anisotropic 



(2+l)D classical XY model 

Z= fvOexp J K {1>v) cos(6 v - 6 X ) 1 , (2.23) 

l(i,i'> J 

where the sum is now over all near-neighbor bonds in both the space and time directions, 
i.e. 1 = (x,y,r). For spatial bonds, 

K l>v = K x , (2.24) 
given by Eq. ( [2.15| ), while for temporal bonds 

K\v = K T , (2.25) 



given by Eq. (|2.22 ). It is implicitly assumed here that the difference between the Villain 



action and the cosine term (which is small for small Ar) is in fact irrelevant in the renor- 
malization group sense. 

Note that we need to take the limit Ar — > which implies K x and K T — > 00 such 
that the geometric mean, 

K = (K X K T ) 1/2 = 1 , (2.26) 



is finite. Universality properties are unaffectedEI if we rescale space and time so that we 
obtain finally an isotropic (2+l)D XY model 



I (MO J 



(2.27) 
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We are interested in the behavior of the boson Hubbard model at T = 0, which means 
taking the number of time slices, M, to infinity. The coupling constant, K, then controls 
the quantum rather than thermal fluctuationsSO. 

Allowing for a non-integer boson density and/or including the random potential in the 
boson Hubbard model, Eq. ( |2.2j ), makes the model more realistic but complicates the sit- 
uation by breaking the particle-hole symmetry of the bosons. This corresponds to broken 
time-reversal symmetry for the quantum rotors (since particle number is represented by 
angular momentum) and hence leads to complex weights in the corresponding classical sta- 
tistical mechanical problem^ which is no longer in the universality class of the (2+l)D XY 
model. The difficulty of complex weights can be avoided by considering the alternative ap- 
proach to Eq. ( |2.19| ) in which we integrate out the {6} variables. Let us do this first for the 
case of integer boson density and no disorder. Adding the effects of disorder and non-integer 
density will then be easy and will lead to a real action. 

We first reexpress the cosine in Eq. (|2.19j) as the best Villain approximation to it, i.e. 



exp(K x cosfl) — ► ex P) ~ 2nm ) 2 \ • ( 2 - 28 ) 



To determine K x we require that the range of the functions on the two sides of Eq. fl2.28| ) 



(as 6 varies from to tt) are the same (the precise angular dependence of the two sides 
will be different but this is presumably irrelevant). Using the Poisson summation formula, 
Eq. ( p. 20 ) and noting that K x — > from Eq. ( |2.15| ), one findsH 



^ = ^ln(|-) • (2.29) 



Inserting Eq. ( |2.28j ) into Eq. (|2 .19 ) and Fourier transforming^ one can carry out the {9} 



integrations exactly. Their effect is to enforce conservation of integer- valued currents defined 
by 

J = (J x , J y , J T ) . (2.30) 

In other words, the current should be divergenceless at every site in space and time, i.e. it 
should obey a continuity equation 

d u J u = . (2.31) 

If J? x , v ,t) ^ es on the bond between sites (x,y,r) and {x + 1, y, r) then it is convenient to 
define J[ x x yT ) = —J(x-i,y,r), etc. The divergence constraint is then imposed at each site by 
requiring that J2v J(t,t) = 0> "where v runs over ±x, ±y, ±r. We thus obtain 

Z^E'expj-^E E M4,r)) 2 } . (2-32) 

{J} [ Z (r,r) v=x,y,r J 

where the sum is over all integer values of the J u from — oo to oo, the prime indicates the 
constraint that J be everywhere divergenceless, and the couplings are 



10 



K T = UAt 



(2.33) 



and K y = K x given by Eq. (|2.29|) . Note that in taking the quantum limit, Ar — > 0, the 
spatial couplings K x and K y diverge, while the coupling in the time direction, K T , tends to 
zero. This is the opposite of what we found in the phase representation, see Eqs. ( 2.15Q and 

We interpret J u as the "relativistic" 3- vector current with ( J x , J y ) being the spatial cur- 
rent and J T being the particle density. Consider the divergenceless current configuration 
represented by the closed loop in the x-r plane shown in Fig. |2|. The physical interpretation 
of this is that at time T\ a boson hops from position x\ to position X2 creating an instan- 
taneous burst of spatial current. This represents a tunneling event in which we assume the 
barrier is high enough that the tunneling time is small compared to the separation between 
time slices in our lattice and hence the event can be treated as instantaneous. This approx- 
imation affects the ultraviolet details of the calculation but is irrelevant to the universal 
zero-frequency behavior. The two vertical lines represent time-like components of the cur- 
rent indicating that there is now a missing boson at X\ and an excess of one boson at x<i- 
After some additional random motion, a boson hops back to the original site, leaving the 
system in the vacuum state at time t^. 

This interpretation of the current is confirmed by consideration of the effect of an external 
vector potential which modifies the potential energy of the quantum rotors to 



V 



(r,r) 



E 

=x,y 



COS 



(2.34) 



where stands for the line integral of the vector potential along the link from site r to its 
neighbor in the t^-th direction. Making this substitution modifies Eq. ( |2.32[ ) with the result 



Z « E'exp \ - £ 

{J} I (r,r) 



E T « +l E J(r, T) A 



(r,r) 



,u=x,y,r 



v=x,y 



(2.35) 



We note that 



J (r,r) 



JlnZ 



v = x,y 



^2.36) 



which means that must thus be the full, physical, gauge-invariant current, not simply 
the paramagnetic piece of the current. 

From our interpretation that J T is the particle density it is now straightforward to include 
both disorder and a value of the chemical potential which gives a non-integer density, and 
one finds 



(r,r) lv=x,y,r 



(2.37) 



We now assume, as in Eq. ( |2.23|) , that the universality class is unchanged if we make the 
couplings isotropic, i.e. 
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J 



exp 



1 

K 



E 

(r,r) 



2 J (r,r) 



s 



^2.38) 



where 



Vr= u 



[2.39) 



and K is a dimensionless coupling constant which has to be adjusted to bring the system to 
the critical point. Varying K corresponds to changing the ratio t/U in the boson Hubbard 
model, Eq. ( |2.2|) , keeping \x/U and A/U fixed. Noting the invariance of the action under 



J T + l 



[2A0) 



we take for simplicity the "largest possible" disorder by choosing \x= 1/2, A = A/U = 1/2. 
The average particle density is then 1/2. Note that this choice of parameters has a statistical 
particle-hole symmetry!! since, upon ensemble averaging the Hamiltonian is invariant under 
the transformation J T — ► — J T , although the presence of the random potential destroys 
microscopic particle-hole symmetry. We have argued above that lack of microscopic particle- 
hole symmetry^ changes the universality class from that of the (2+l)D XY model, so one 
can ask whether having statistical particle-hole symmetry changes the universality class from 
that of the generic Bose glass to superfluid transition. At least in one dimension, the answer 
is no, as shown by FisheiH, and we shall assume that the same is true in d = 2. 

We have thus arrived at a representation of the original quantum problem, involving 
integer link variables. Noting that the J T represent the boson density, it can be thought of 
as an imaginary time "world-line" path-integral representation of the problem, simplified to 
the extent that it treats just the phase fluctuations of the underlying Hamiltonian. 

The partition function, Eq. (|2.38|) can be written in terms of an effective (2+l)D classical 
Hamiltonian or action, given by 



H 



v 



1 

~K 



E 

(r,r) 



r (r,r) - (V + Vr)J{r, T ) 



(2.41) 



Evidently, when v r = 0, integer values of /2 can be absorbed into the definition of JL T ), so 
the model reduces to the (2+l)D Villain model, which is in the same universality class as 
the (2+l)D XY model. These points, are, however, just special multicritical points and the 
generic behavior is not that of the XY mo defl. 

We have already noted that the time component, JJ TT \i of the link variables corresponds 
to the particle density or boson occupation number. Long-range Coulomb forces can then 
be introduced in the following way. 



H = Hy + Hq 

e* 2 

ff c = yEE ( J (r,r) - n )G(r - r')(J[ r , iT) - n ) 

T <r,r'} 



(2.42) 
(2.43) 
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Here e* is the effective boson charge, and no, which represents the compensating background 
charge, is the average particle density, and G is the Coulomb interaction. In our simulations 
with long-range interactions, the particle number was always kept constant, as opposed 
to the case where only short-range interactions were present where we always allowed the 
particle number to fluctuate. Calculations of the Coulomb interaction, G(r) must allow for 
the finite lattice size and periodic boundary conditions. We do this by the usual Ewald 
methodil'0. Another way is by a lattice Green's function: 

r(r) = — \^ cos ( k • r ) (o AA) 

{ 1 L 2 ^ [4 — 2 cos(^) - 2 cosfo,)] V* ' ^ ' > 

where k = (27r/ L)(n x ,n y ), with n x ,n y = 0, ... ,L — 1, The term with k = is removed 
to ensure charge neutrality. For large distances and large lattices the Ewald sum and the 
lattice Green's function become almost identical and approach 1/r. However, close to the 
origin the two forms are somewhat different. If the critical properties are universal they 
should not depend on the specific form of the potential close to the origin. We use this as a 
test of our computer codes and of the universality of our results. Indeed as we shall see the 
two forms yield equivalent results. 



III. SCALING THEORY 

In order to better understand the universal features of the phase transition it is very 
useful to consider the scaling behavior of various physical quantities in the regime of the 
diverging correlation length. Such considerations not only tell us why the conductivity is 
universal but will tell us how to analyze experimental and Monte Carlo data to determine 
that one is actually in the critical (scaling) regime. 

From now on, we shall denote the number of time slices by L T , rather than M, so the 
space-time lattice is of size L x L x L T . Periodic boundary conditions will be applied. Note 
that the ground state energy density of the original 2D quantum problem is related to the 
free energy density of the (2+l)D equivalent classical problem since 

/ = -&«(^ lnZ = -F 1 " Tte "- < 3J) 



where H is given by Eq. ( |2.42|) , V = (aL) 2 L T Arh is the "volume" of the (2+1) D space-time 



system, with a the lattice spacing in the spatial directions and hAr the lattice spacing in 
the (imaginary) time direction. 

Since space and time are not equivalent we have two correlation "lengths", £ in the 
space direction and £ T in the time direction. These two correlation lengths will diverge 
with different exponents at the critical point and we can define the dynamical exponent, z, 
through the relation 

£~r", £ T ~£*, (3.2) 

where 6 measures the distance from the critical point, K c , i.e. 5 = (K — K c )/K c . There 
is a microscopic frequency, u c , related to the lattice spacing hAr, in the time direction by 
u> c = 27i/(hAr), so we can relate £ T more precisely to £ as 
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(3.3) 

where b is a microscopic length of order the lattice spacing, a. 

A. Stiffness 

First we discuss the scaling theory describing the singular behavior of the free energy 
density near the critical point. From Eq. ( |3.1| ) one sees that f/H has dimensions of inverse 
(length^ x time). Hyperscalingjlll states that multiplying the singular part of this free energy 
density, f s , by the (2+l)D correlation volume, £ d £ r , one obtains a constant, A say, of order 
unity as the critical point is approached, i.e. 

jt% = A . (3.4) 

In this section, we will frequently give results for arbitrary space dimension, d, even though 
we are ultimately interested in the case of d — 2. One can consider A to be a critical ampli- 
tude for a dimensionless quantity (or combination of quantities) which is finite at criticality. 
According to two-scale factor universalityE^TL^, such quantities are not only constants, but 
are also universal. 

We next discuss the scaling of the extra free energy cost to impose a twist on the phase of 
the condensate. We will use this to locate the critical point, K c , to high accuracy. The extra 
free energy density is related to the superfluid stiffness, also called the helicity mo dulusS, 
which is proportional to the superfluid density of the system. A uniform twist in the phase 
of the order parameter can be introduced by applying a twist of size 6 at the boundary, in 
(say) the x direction. This will then give rise to a phase gradient 

V0 = G/(aL) . (3.5) 



8f. 1 _,™ N 2 



The (zero-frequency) stiffness, p, is then defined bj 

-oiyef , (3.6) 



h 2' 



so 



P = ^™- (3-7) 



h <96 2 

Since G is dimensionless, p has dimensions of inverse (length^ -2 x time). Hence using 
hyperscaling and two-scale factor universality we obtain 

p£ d ~% = C , (3.8) 
where C is another universal constant. Consequently, 

p ~ £-(*f*-2) ; (3.9) 

which is a generalization of the Josephson scaling relation for the classical transition, p s = 
(m/h) 2 hp ~ £ d ~ 2 . The difference is that d is replaced by d + z for the quantum transition. 
This replacement also holds for other hyperscaling relations (i.e. those scaling relations 
involving the space dimensionality). 
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B. Conductivity 



We can extend the notion of a superfluid stiffness, p, to a frequency- dependent stiffness, 
p(iuj n ), where u n = 2imk,BT /% is the Matsubara frequency. The conductivity is then related 
to p{iuo n ) by the Kubo formula^ 

a{iu n ) = 2kG q P -^ , (3.10) 

where Gq = Rq 1 , with Rq defined in Eq. ( pLlD , is the quantum of conductance. The 
quantity p defined in the previous section is given by p = p(0). We emphasize that p is the 
stiffness and not the resistivity. Close to the critical point we can generalize Eq. ( ft.8|) to 
finite frequency^ by the following scaling assumption 

P (tCO n )=e- d C 1 p(^r) ■ (3.11) 

Since the argument of the scaling function p is dimensionless, and so has no non-universal 
metric factors associated with it, the entire scaling function p(x) is universal. Clearly, 
p(0) = C, the same universal constant that appears in Eq. (|3.8|) . Furthermore, since p{iuj n ) is 
finite at finite frequency even at the critical point, one must have, for large x, the asymptotic 
behavior 

p(x) = Dx id+z - 2)/z , (3.12) 

where D is again universal, in order that the dependence on £ and £ T cancels at criticality. 
Substituting this into Eq. ( |3.10| ) and noting Eq. (|3.3p, one has, at criticality, 

Oj \ (d-2)/z 



({. \ (a— z)/z 
— . (3.13) 

UJ r .J 



Immediately we see that when d — 2 all microscopic lengths and frequencies drop out so the 
d.c. conductivity is universali at the critical point, given only by fundamental constants and 
the universal dimensionless number D. The universality of the d.c. conductivity is analogous 
to the universal jump@ in (h/m) 2 p s jksT c at the finite-temperature Kosterlitz transition!. 
In fact this quantity corresponds, essentially, to Eq. (|3.10| ) with hu n replaced by ksT c . 
Strictly speaking Eq. ( 3.13|) only refers to the singular part of the conductivity. However 



since we approach an insulating phase where the conductivity must be zero, the conductivity 
cannot have an analytic part at the critical point. 

C. Compressibility 

The compressibility, k, is defined by 
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where n is the boson density and fi the chemical potential. We can write an expression 
equivalent to Eq. ( |3.6|) for the compressibility by notingS that the Josephson relation (for 
imaginary time) is 5[i = 86 /dr, so 

5f= 1 -K(d T 6)\ (3.15) 

i.e. we apply a twist in the (imaginary) time direction instead of along one of the space 
directions. Note that it is the total compressibility which enters this expression!!!. Following 
the arguments that led to Eqs. (|3.8|) and ( |3.9p one finds 

^V = ~£V0V = const. (3.16) 

so 

£-(«*-*). (3.17) 

Fisher et al.0 have argued that z — d at the Bose glass to superfluid transition and so the 
compressibility is finite at criticality. We shall see that our numerical results support this. 
Note that even in this case, the compressibility is non-universal at criticality, because the 
non- universal factors b and uo c appear in Eq. (|3.16|) . One can also determine the form of the 
wave- vector dependent compressibility at criticality, following the scaling arguments that we 
used above to determine the conductivity. One finds 

«(*;) ~ k d ~ z . (3.18) 



IV. QUANTITIES OF INTEREST AND FINITE SIZE SCALING 

In this section we show how to calculate the quantities of interest from the Monte Carlo 
simulations, and we discuss the finite size scaling techniques which we will need. Having 
demonstrated in the last section, that the lattice spacings, a and HAt, do not enter expres- 
sions for universal quantities, such as the conductivity at the critical point, we set these 
lattice spacings (and ft) to unity from now on. 

To perform the quenched disorder averages it is necessary to first do a "thermal" average 
over the J u variables, denoted by (•••), for a fixed realization of the quenched disorder 
potential, and then average observables over the quenched disorder v r , which we indicate by 



A. Stiffness 

To calculate the uniform stiffness, p(0), note from Eq. ( |2.34| ) that a uniform twist in 
the x direction becomes equivalent to considering the system in the presence of an external 
vector potential of the form 

A x r = d x B8 x , v . (4.1) 
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From Eqs. flOj ), Q and (0) one finds that 



P(0) 



L d L T 




(4.2) 



Near the critical point, the correlation length is much larger than the size of the system, 
so finite-size effects will be important. We therefore need to derive a finite-size scaling formH 
for the stiffness. The basic finite-size scaling hypothesis is that the size of the system only 
appears in the ratio L/£, and, for quantum problems, the corresponding ratio in the time 
direction, L T /^ T . Thus we have 



p(o) = e 



-(d+z-2) 



which can be more conveniently expressed as 

1 



P(0) 



L d+ 



2-2 



(4.3) 



(4.4) 



where P and p are scaling functions. It is thus essential to work with system shapes for 
which the aspect ratio 



(4.5) 



is a constant, otherwise the scaling function p depends on two variables and is complicated to 
analyze. If this is done, L d+Z ~ 2 p is independent of L at the critical point 5 = 0. Furthermore, 
in the disordered state, the system is insensitive to changes in the boundary conditions if 
the size is bigger than the correlation length, so L d+Z ~ 2 p will decrease (exponentially) with 
increasing L. By contrast, in the ordered state, p tends to a constant so L d+Z ~ 2 p increases 



with increasing L. Thus, the critical point is located at the intersection of curves for L d+Z ~ 2 p 
as a function of coupling K for different lattice sizes. One can then determine v from Eq. ( |4.4|) 
by requiring that the data for different sizes (but fixed aspect ratio) collapse on top of each 
other in a plot of p(0)L d+z ~ 2 against L l / U 5. Note that in order to choose the sample shapes 
in the simulation, we need (unfortunately) to have already made a choice for z. 

Since the current is divergenceless, we can divide the configurations into different topo- 
logical classes according to the winding number of the boson world lines around the torus 
of size L in the space direction 



n. 



J (r,r) > 

(r,r) 



(4.6) 



so that the stiffness is simply proportional to the mean-square winding number 

P(0) 



1 



n, 



(4.7) 



It is instructive to comment on the analogy to the Feynman ring exchange picture of super- 
fluidity in liquid helium^. Rather than viewing a nonzero winding of a boson world line as 
an event involving a single boson, we can view it as formed by adding up a chain of hops 
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of many bosons. This is closely analogous to a Feynman ring exchange, and adds some per- 
spective on the transition in our model. The superfluid stiffness arises due to macroscopic 
condensation of ring exchanges of global currents carrying nonzero winding number, and 
the critical point is where the gain in free energy from the "entropy" of the ring exchanges 
matches their energy cost. In the presence of macroscopic ring exchanges which wind around 
the sample, the free energy is sensitive to Aharonov-Bohm flux (boundary condition twists) 
and hence the system exhibits off-diagonal long-range order in the conjugate phase variable. 



B. Conductivity 



The frequency-dependent stiffness involves the Fourier transform of the current-current 
correlation function 

pl*>n) = Jj^Kl E ^" T Jf rjr )i 2 }]av , (4-8) 
T (r,T) 

where, with the lattice spacings in the space and time directions set to unity, the Matsubara 
frequency is given by u n = limj L T , and r is now an integer, 1 < r < L T , denoting a 
particular time slice. In these units, the conductivity is still given by Eq. ( p.lOj ). 



C. Compressibility 

From Eq. (|3.14 ) it follows that the zero wave vector compressibility is given by 

«(0) = J^[(K> ~ WV , (4.9) 
where Nj, is the total number of particles, 

^ = J" E J (r,r) ■ (4-10) 
(r,r) 

The last term in Eq. ( |4.9Q involves the square of a thermal average. This term is thus likely 
to give systematic error® if determined within one replica, so we evaluate it as [{N a )(Np)} av , 
where the indices refer to two different replicas. 

If global moves are not included, the boson density is a constant, and consequently k, as 
defined, is zero. However, the wave-vector dependent compressibility n(k) is nonzero, and 
one can obtainlll, estimates of k = k(0) by taking the limit /c — *> 0, even when global moves 
are not performed. 

The finite-size scaling form for the compressibility follows from arguments similar to 
those used above for the stiffness and is 

' r LV»6, . (4.11) 



K = ; K 



L d ~ z V ' L 1 
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D. Correlation functions 



Consider the following correlation function 

C(v,v',T,r') = [(e^-^^% v , (4.12) 

where the 9's are operators for for the phase of the bosons, and e l9r( - T ' 1 = e TH e l9r e~ TH . We 
shall see that this correlation function gives information on a third critical exponent, 77, 
defined in Eq. ( 4.18Q below, in addition to the exponents v and z already discussed. Due to 



translational invariance, C(r, r', r, r') = C(r — r', r — r' 

We shall here only consider two basic types of correlations. The equal-time correlation 
function where r = 0, r = (x, 0), and the time- dependent correlation function at time r 
and r = 0. By redoing the argument which led from Eqs. ( |2.4| ) to ( |2.32| ) for the correlation 
function rather than for the partition function, one finds that the equal-time correlation 
function can be expressed as 

C x (r) = [ ( I] exp {~{Jf r ,r) + h) > lav , (4-13) 

r Gpath 

where "path" is any path on the lattice at fixed r connecting two points a distance r apart 
along the x-direction. For each link on the path v = x or y, depending on whether the link 
is along the x or y direction. The simplest case, which was used in the simulations, is the 
straight line path, in which case all the link variables in Eq. ( |4.13| ) are J x . A very similar 
result is is found for the usual Villain modelil. Since the Hamiltonian is invariant under a 
sign change of J x , it follows that C x {r) = C x (—r). 

We now turn to the time-dependent correlation function, 

C+(t) = C(r = 0, r) = [<e^M-^°»)] av . (4.14) 

Physically this is the Green's function for creating a particle at imaginary time and de- 
stroying it at time r. This correlation function can be expressed in terms of the link variables 
in the following form 



C + (r)=[( J] e X p\--(- + Jl rtT) -fi t 

t Gpath 



K 



EWrV) " n )G(r - r') + ^ £(G(0) - G(r)) 



(4.15) 



In this expression, "path" is the straight line path between two points with the same space 
coordinate, r, starting at imaginary time equal to 0, say, and ending at a later time r. A 
more general expression for a path wandering in the space directions can also be derived. 
One can also consider 

C_(r) = [<e-^-^°»}] av , (4.16) 

which is the Green's function for creating a hole at imaginary time and destroying it at r. 
In terms of the link variables 
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T 



[( n exp 

r gpath 

e* 2 r ' 



-(- 



J (r,r) 



n )G(r 



2L 2 



E(G(0)-G(r)) 



(4.17) 



Except when there is statistical particle-hole symmetryEl, C + (r) 7^ C_(r)H. However, 
one can show that C+(t) = C-(L T — r), which corresponds to the equivalence between a 
particle traveling forwards in time and a hole traveling backwards. This will be useful in the 
simulation because the statistics get worse with increasing r so, for r > L T /2, it is better 
to compute C_(L r — r) than C + (r). As required, the correlation functions are periodic, i.e. 
C(t) = C(t + L T ) for both C_ and C+. 

Following Ref. [Il] we now make the assumption that the long-distance, large-time be- 
havior of the correlation functions will be given by the scaling form 



C(r,r) = r-^- 2+ ")/(r/e,r/r) 



(4.18) 



which defines the exponent rj. If r approaches zero but r remains finite, the correlation 
functions must remain finite and nonzero. Thus we obtain 



C(r = 0,r) = r -( rf +^- 2 +")/^( r /f 



(4.19) 



At the critical point we should therefore have 

CJr) ~ r 



•Vx 



CJt) 



T 



-Vt 



(4.20) 



where 



y x = d + z- 2 + r] 

y T = (d + z-2 + 7])/z . 



(4.21) 



Thus the power law fall-off of the correlation functions at criticality determines both 77 and 
z provided the correlation functions can be evaluated for large enough system sizes that 
finite-size corrections are unimportant. 



V. MONTE CARLO METHODS 

To satisfy the zero divergence criterion in Eq. ( |2.31| ) our basic (local) Monte Carlo move 
consists of changing all the link variables around one plaquette simultaneously in the manner 
shown in the lower left corner of the Fig. |3], thus changing the local current. Two of the 
link variables are increased by one, the other two decreased by one. An equivalent move 
going in the other direction is also used, i.e. the plusses and minuses are interchanged. In 
addition, we need to include non-local moves to fully equilibrate the system. The global 
moves consist of changing by ±1 a line of link variables stretching all through the system. 
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Nonlocal moves are included in all three directions 5 = x,y, r, except when the model has 
long-range interactions, in which case no global moves in the time direction are performed 
in order to keep the particle number constant. Global moves in the time direction amount 
to either introducing or destroying a boson. It is easy to see that global moves in the 
space directions correspond to a change in the winding number^E, defined in Eq. fl4.6| ). 
The nonlocal moves we use are illustrated in Fig. |3|. One Monte Carlo sweep of the lattice 
consists of a sweep of local moves followed by a sweep of global moves. 

Due to the continuity equation, Eq. ( |2.31| ), the sum of at a given time slice, J2r ^ri 
is always the same for any value of r, albeit this constant may vary as a function of Monte 
Carlo time because of global moves in the time direction. Likewise, the sum of J x in any 
y-r plane will be the same for all such planes at a fixed Monte Carlo time, and similarly for 
the sum of the J y . 

Expectation values of observables have to be computed by quenched disorder averaging, 
which is known from the study of spin glasses^ to have many potential pitfalls. Close to 
the critical point we typically have to average over from 200 to 1000 different realizations 
of the disorder, and somewhat fewer away from the critical point. It is crucial to carefully 
assure that the J u variables are thermally equilibrated. The equilibration time at the critical 
point for our update scheme varies with system size L as r mc ~ L Zmc , where z mc is the Monte 
Carlo dynamic exponent. For the short-range interaction case we have determined^! z mc ~ 6 
so that extreme caution is required in attempting to equilibrate large lattices. We take an 
approach similar to what has been done for spin glass systemsEa Two identical replicas are 
run in parallel for a given realization of the disorder. We define the "Hamming" distance 
between replicas a and (3 as: 

KM = E [ J {r,T),a(t0 +t)~ 4, T ,^(t0 + t)] 2 , (5.1) 
(r,r) 

where to is the number of Monte Carlo sweeps (MCS) used for equilibration, and t is the 
number of subsequent MCS. We also define a "Hamming" distance for one replica at two 
different Monte Carlo times, 

Kit) = [J(r,r),a(t + k) ~ J(r,T), a (t )] 2 . (5.2) 

(r,r) 

We determine \h v a n(io)]av an d [^o(^o)]av for a sequence of values of to increasing exponen- 
tially, t = 10, 30, 100, 300, 1000, . . ., up to t = T , so T is both the number of MCS for 
measurement and the number of MCS for equilibration. If t is sufficiently large that the 
system has equilibrated, one has [Kai^o)]™ = IAa(*o)]av> an d we made sure that this con- 
dition was fulfilled, at least for to = To. To achieve equilibration we took To to be of order 
3,000 for the smaller system sizes but found that we needed up to 30,000 for the larger sizes. 
Since the different disorder realizations give statistically independent thermal averages, we 
can estimate the statistical error from the standard deviation of the results for different 
samples. Note that there are big sample to sample fluctuations, so it is necessary to average 
over a large number of samples. In order to study as many samples as possible within the 
available computer time, we only run each sample for the minimum number of MCS neces- 
sary to get a few statistically independent measurements. This is why the number of sweeps 
for averaging is the same as the number used for equilibration. 
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VI. SHORT RANGE INTERACTIONS 



In this section we shall assume that no long-range Coulomb interactions are present. 
Furthermore, we shall always take the random chemical potential to be specified by 

A = i A = i. (6.1) 

The reason for this choice of jl is that we want to be as far away as possible from any Mott 
insulator phased, and these are centered on integer values of jl for weak disorder. The choice 
of A was influenced by the need to make the disorder not too small (otherwise the effects of 
disorder would only be seen for large sizes which we are unable to simulate) and also not too 
large, because this effectively makes U small, and so, again, the asymptotic behavior may 
only set in for large sizes. In the absence of more detailed information, it seems sensible to 
make all the important couplings of comparable size. We again emphasize that universal 
quantities like the critical conductivity are independent of these details. 

Some inequalities involving the critical exponents, u, rj, and z have been obtained. First 
of all, Fisher et al.0 have argued that the compressibility is finite at the transition and so 

z = d . (6.2) 

Fisher et al.0 also argue that 

rj<2-d, (6.3) 

on the grounds that the density of states should diverge as the transition is approached 
from the Bose glass side. In addition, since the correlations must decay with distance at 
criticality, it follows from Eq. fl4.18|) that 



d + z-2 + T]>0. (6.4) 

Note that since z > one can have a negative _n even in two dimensions. There is also a 
general inequality applicable to random systems§lll 

v > 2 - , (6.5) 

which is a generalization of the Harris criterion^. The value of the dimension that should 
be inserted into this expression is the number of dimensions in which the system is random, 
i.e. the space dimension d and not d + 1 or d + z. 

As noted in the discussion below Eq. ( [4.5[ ) we need to know the dynamical exponent z 
in order to choose sample shapes which allow a simple finite-size scaling analysis, i.e. the 
samples should be of size L x L x cL z , where c is the aspect ratio. Most of the simulations 
were done assuming z — 2, the value predicted by Fisher et al.0. We have done additional 
simulations with shapes corresponding to other values of z, but find that the scaling is much 
less good if z is significantly different from 2. For z = 2 we have taken two different aspect 
ratios 1/2 and 1/4, with the following systems sizes: 4x4x8, 6 x 6 x 18, and 8 x 8 x 32 
for aspect ratio 1/2, and 6x6x9, 8 x 8 x 16, and 10 x 10 x 25 for aspect ratio 1/4. We 
were unable to study larger lattices because the relaxation times were too long. 



As a test of our program we checked that we were able to reproduce the results of Ref. ^0 
in the absence of disorder and with jl = 0. We found complete agreement!^ between the two 
simulations. 



22 



A. Equilibration 



We test for equilibration using the method described in Section M. As an example, Fig. |4] 
shows the Hamming distance for the x and r link variables for a system of size 8 x 8 x 16, at 
the critical point. We see that the system equilibrates rather quickly in about 1,000 MCS. 
Also, we see that the x and r link variables equilibrate in roughly the same time, as one 
would expect since they are coupled through local moves. 



B. Determination of the Critical Point 

We start the analysis by locating the critical point. Since, as discussed above, we assume 
that z = d = 2, the relevant quantity to plot, according to the finite-size scaling analysis 
in subsection [IV A| , is p(0)L 2 . Results for aspect ratio 1/4 are shown in Fig. |5|. Since 



the critical point is located where the curves cross, the figure demonstrates clearly that 
there is a transition close to K = 0.25 between a superfluid phase for K > K c with finite 
superfluid density, p s (remember that p s ~ p(0)), and an insulating phaseafor K < K c with 
zero superfluid density. Our best estimate of the critical coupling is K c = 0.248 ± 0.002. 
A substantial amount of computation went into the production of this figure. Close to the 
critical point 1000 to 2000 disorder realizations were performed, with, for the largest size, an 
equilibration time of T = 10, 000 followed by 10,000 MCS for averaging with a measurement 
every 10 MCS. 

Simulations with aspect ratio 1/2 were also performed and the same critical coupling 
was found, as expected since this is a bulk property. 



C. The Compressibility 

We now turn to the compressibility. Fig. |^ shows the compressibility, as calculated from 
Eq. (|4.9|), for a range of different couplings centered around the critical coupling K c = 0.248, 
for lattices with aspect ratio 1/4. We see that the compressibility remains finite through the 
transition, including in the insulating phase, K < K c . This is consistent with the prediction 
that the insulating phase should be a Bose glass with finite compressibility in the presence of 
disorder^. According to the scaling theory, Eq. (|3.17|) a finite compressibility at criticality 
implies z — d (— 2), as argued by Fisher et al.M. By contrast, simulations performed^ with 
no disorder and p, — 0, where the model becomes equivalent to a (2+l)D XY model, find 
that the compressibility vanishes in the insulating phase, consistent with it being a Mott 
insulator. 

Fig. shows the wave- vector dependent compressibility for the aspect ratio 1/4, at the 
critical point K c = 0.248. Similar results have been obtained for the aspect ratio 1/2. 
Clearly there is no dependence on the wave vector as expected from Eq. ( |3.18| ) and the 
result z = d. 

From the above we have established that the insulating phase is indeed a Bose glass and 
not a Mott insulator, at least for the strength of the disorder that we have been considering 
here, A = 1/2. This is in agreement with previous studiesSSSil. 
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As further evidence of the existence of the Bose glass we now turn to a discussion of 
the correlation functions in the insulating phase. One important prediction of the scaling 
theoryEl is that the Green's function should vary with a power of imaginary time, 

C(t = 0,t)~ P i(P)/t, (6.6) 

rather than exponentially, as might have been expected. Here pi(0) is the single-particle 
density of states at zero energy. In order to check this prediction we did simulations deep 
in the insulating phase, K < K c . Fig. |8] shows the time- dependent correlation function, 
C + {t) = C(r = 0,r), for a system of size 8 x 8 x 16 at a coupling equal to K — 0.175, 
well below K c ~ 0.248. The right hand part of the figure is obtained by calculating C_(r), 
and using the relation C + (r) = C_(L T — r) discussed in subsection [IV D| . For each disorder 
realization, 30,000 MCS were performed, followed by another 30,000 MCS to do the thermal 
averaging, and finally we averaged over 100 different disorder realizations. The relatively 
elaborate thermal averaging was done in order to obtain small error bars at large r. The 
dashed line is a power-law fit to the form 0. 170(2) (r~ L1 °( 8 ) + (L T — r) -1 ' 10 ®) where the 
numbers in parentheses indicates uncertainties on the last digit. This fit used all data points 
shown and gave a goodness of fit of 0.84 and \ 2 — 7.9. Here we define the goodness of fit to 
be T((N — 2)/2, x 2 /2), where N is the number of data points and T is the incomplete gamma 
function. No sign of an exponential dependence on r was observed. The errors indicated are 
statistical and do not include possible systematic errors. A fit at K — 0.15 yielded a similar 
value for the exponent of 1.05 ± 0.04. We conclude that the time-dependent correlation 
functions clearly display power-law behavior in the Bose glass phase and, furthermore, the 
associated exponent is close to 1 as predicted by scaling theory^. 



D. The Conductivity 



In the thermodynamic limit, L — > oo, at vanishingly small T (L T — > oo), and for u n — > 0, 
the conductivity at the critical point should tend to a finite, universal value, a*, as discussed 
in subsection [111B| . In the simulation there will be various corrections to this. First of all, 
one might ask whether the order of limits T — > and u n — > affects the value of a*, even in 
the thermodynamic limit. For the case of no disorder and integer filling, where the transition 
is to the Mott insulator, the answer is certainly yesS. In this case the conductance is finite 
if the T — > limit is taken first whereas a* = oo if one first takes the zero frequency (d.c.) 
limit because a persistent current can flow in the absence of umklapp processes, which vanish 
as T — > 0. However, in the presence of disorder, the d.c. conductivity is finite as T — > and 
so we see no reason why the order of the limits, T —>■ and u n — > should play a role for 
transition to the Bose glass phase discussed here. We shall see that allowing for a dependence 
on u) n /T oc LU n L T does give a slighly better fit for the case of short-range interactions, but the 
value of a* is not changed significantly. Given the rather limited range of sizes that we can 
study, we feel that the results are consistent with there being no dependence on u n /T. One 
might also be concerned about finite-size corrections to the conductivity, but our results are 
consistent with their being very small. Of course, the conductivity is frequency dependent 
and so will differ from the universal value when u n becomes comparable with some other 
scale, such as the ultraviolet cutoff set by the lattice spacing. 
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Fig. |9] shows the resistance per square (which is the same as the resistivity) plotted against 
frequencyH, evaluated from Eq. (p. 10 ) for aspect ratio 1/4 at the critical point K c = 0.248. 



Again considerable computation has gone into the production of this graph in order to obtain 
good statistics. For the two smallest system sizes about 2,000 disorder configurations were 
generated while for the largest size only 1,000 were done. From 3,000 to 10,000 MCS were 
done for equilibration followed by the same number of sweeps for measurements. The data 
collapse is excellent. 

To determine the universal conductivity we have to analytically continue the MC data to 
real frequency and extrapolate to uj = 0. For typical quantum MC simulations, this analytic 
continuation is extremely difficult to perform. However it turns out to be straightforward in 
the present case, since the data for the resistivity varies linearly at small u n , which implies, 

a* 

= 77—1 — • ( 6 - 7 ) 

(1 + \u n \T ) 

This is easily seen to analytically continue to the Drude form of the conductivity: 

<T(u + i6) = T a * . (6.8) 
(1 - iivtq) 

Thus the boson system at the critical point is neither insulator nor superfluid but rather a 
Drude metal. The Drude parameter tq ~ 1/uj c is a non-universal relaxation time controlled 
in our model by the ultraviolet cutoff. 

Assuming this linear variation of the resistivity with u n , a least squares fit has a very 
small error. The main source of error in the determination of the d.c. conductivity therefore 
comes from the uncertainty in the determination of the critical point. We estimated this 
error by making the same linear fit to the resistivity data at the ends of the interval given 
by the error bars of the critical coupling. From all our data for two different aspect ratios 
we finally estimate 

c* = (0.14 i 0.03)Gq , R* D = (7 A ±1.6) R Q . (6.9) 

The universal conductivity has previously been calculated for the case of no disorder 
and integer boson fillingll, where the insulating phase is a Mott insulator, rather than the 
Bose glass discussed here. In that case a* ~ 0.285Gq. Thus we find that even though the 
present model is somewhat more realistic, including the disorder takes us further from the 
experimental value which is in the vicinity of unity. A suitably defined universal conductance 
can be calculated exactly in 1130 both with and without disorder. The exact solution in ID 
shows that the ratio of a* in the dirty case and in the pure case is exactly 3/4. This is of 
the same order of magnitude as the ratio 0.14/0.285 ~ 0.5 between the MC results in 2D. 
Hence we see that the trend of decreasing critical conductivity upon adding disorder is the 
same as for the exact solution in ID. 



E. The Exponent v 



To determine the correlation length exponent we try to collapse the data in a scaling 
plot of p(0)L 2 versus SL 1 ^, based on Eq. (|4.4j). The plot is shown in Fig. ITDL for which the 
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parameters used are K c = 0.248 and v = 0.9. From this and other plots for different values 
of the aspect ratio we estimate 

v = 0.90 ±0.10 . (6.10) 

Interestingly, the inequality, v > 2/d, (with d = 2 here) derived by Chayes et al.il is only 
just satisfied and may, in fact, be an equality for this model. The equality, v = 2/d, has 
been found for certain models of correlated disorderl!. 



F. The Correlation Functions 

So far, we have determined the universal values of v and a*. We have also found that 
the finite-size scaling works best with z = 2 and, according to the scaling theory, our results 
for the compressibility agree with this. In this subsection we conclude our discussion of the 
model with short-range interactions by looking at the correlation functions, which give us 
the value of the third exponent, 7], and another estimate for z. 

The data for C x (r) for sizes L = 8 and 10 with aspect ratio 1/4 at K c = 0.248 are shown 



in Fig. 11. In order to make full use of all the points we fit the correlation functions to the 



following form 

C x (r) = c(r~ y * + (L - r)~ y *) , (6.11) 

which takes the periodic boundary conditions into account. For L = 8 the best fit had the 
form0.18(l)(r- 2 - 02 « + (L-r)- 2 - 02 «). For L = 10 the fit was 0.18(l)(r- L94 ( 2 ) + (L-r)- L94 ( 2 )). 
These fits are indicated by the dashed lines in Fig. 0. 

In order to obtain rj and z from Eq. (f4.21|) we also need to obtain results for the time- 



dependent correlations. Since fi = 0.5 there is statistical particle-hole symmetry^ and so 



C-(t) = C + (t). Fig. [12| shows the results at the critical point K c = 0.248 for an aspect 
ratio of 1/4, where the data for r > L T /2 were obtained from C_(L T — r). The linear sizes 
were L = 6 and 8. We fit to the form 

c(r"^ + (L T - r)-^) , (6.12) 

and we find for L = 6 the optimal fit has the form 0.266(4)(r- L03 W + (L - r)~ L03 W). 
For L = 8 the best fit has the form 0.256(4) (r" - 94 ^) + (L — r)- a94 W). We see that the 
exponent governing the power-law behavior is about 1/2 of the equivalent exponent in the 



space direction, indicating from Eq. (4.21) that the dynamical exponent z must be close to 



2. Combining all our estimates for z we find 

z = 2.0 ±0.1. (6.13) 
Our estimate for r\ obtained from Eq. (f4.21|) , including results from the two aspect ratios, 



is 

7] = -0.1 ±0.15 . (6.14) 

This agrees with the inequality Eq. ( |6.3j ) r] < 2 — d (= 0), which is possibly satisfied as 
an equality. Noting that v is given by Eq. ( |6.10|) and a* by Eq. ( |6.9|) , this concludes our 
discussion of the universal properties of the short-range model. 
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VII. LONG-RANGE COULOMB INTERACTIONS 



We shall now include long-range Coulomb interactions. Throughout this section we again 
take 

A = i A = i, (7.1) 

but, in addition, take the charge of the bosons to be nonzero. Most of our studies used 

e* 2 = ~ , (7.2) 

but we also have some results for e* 2 = 1/4 as a check that the strength of the Coulomb 
term is irrelevant. The long-range interactions force us to keep the total boson number fixed, 
since it has to be compensated by a (fixed) background charge to avoid an energy which is 
infinite in the thermodynamic limit. We do not, therefore, allow global moves in the time 
direction (the r-link variables). 

As in the section above on short-range interactions we first discuss what inequalities and 
estimates there are for the exponents. The result z = d, quoted earlier is only applicable to 
short-range interactions. For a 1/r potential, FisheiS has argued that 

z = l . (7.3) 

A simple way to see this§ is to compute the characteristic energy, Ae, given by the potential 
energy at r = £, i.e. Ae = G(r = £) ~ where G is the 1/r Coulomb potential. If this 
is the relevant energy scale in the problem then As ~ with z = 1. This argument is 
trivially generalized to interactions falling off with some arbitrary power of the distance, 
G{r) ~ r~ A , and leads to z — A. We expect that this is valid for A smaller than d, the value 
for short-range interactions, and that for larger A, the the dynamical exponent sticks at its 
short-range value, z = d. Based on scaling of a renormalized charge Fisher et al J derive the 
inequality 

z < 1 , (7.4) 

and an argument that the second sound velocity should not diverge at the critical pointlll 
gives, quite generally, 

z > 1 . (7.5) 

Hence there is quite strong evidence that z = 1 for the 1/r interaction. This is convenient 
for the Monte Carlo work, because, although the computer time per update increases by 
adding the long-range interaction, the number of lattice points is not so large as for the 
short-range case because we only have to scale L T with the first power of L, rather than its 
square. 

In the section above on short-range range interactions, the inequality rj < 2 — d was 
discussed. This was derived0 with the assumption that the density of single-particle states 
in the Bose glass phase is finite at zero energy, an assumption which is no longer correct 
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with 1/r interactionsBSEE In a classical model, called the "Coulomb glass", Efros and 
Shklovskiffi [ES] argue that that the single-particle density of states, p\{e), vanishes at 
e = 0, due to the "Coulomb gap". Assuming that 

Pi{s)~e\ (7.6) 

ES obtain a bound on the density of states for small £, 

Pi(e) < C£ d - 1 , (7.7) 

so that 

a>d-l. (7.8) 

The value of a in 2D does not seem to be precisely knownB. 

Since the Coulomb glass model is classical, the statistics of the particles does not matter, 
so Eq. (|7.8|) should be applicable to the Bose glass phase, provided quantum fluctuations are 
unimportant in this region, as is argued for the electron caseHS. It is then straightforward 
to determine the long time behavior of the Green's function in the Bose glass phase. For 
t > we have 

CM = r d£e- £lrl Pl (£) (7.9) 
Jo 

1 



7-1+a 



The same argument, together with Eq. (EL21|) , indicates that at criticality pi(e) ~ ^ d - 2 + r i)l z . 
We now assume, following Fisher et al.&3, that pi{£) at small e grows as the critical point 
is approached, in order to match onto the delta-function density of states in the supercon- 
ducting state. In other words, the density of states exponent is smaller at the critical point 
than in the Bose glass phase, i.e. (d — 2 + rj)/z < a, or 

7]<2-d + az , (7.10) 



which, as expected, reduces to Eq. ( p73|) for a constant density of states, a = 0. The bound 
v < 2/d, Eq. ( |6.5| ), should also be valid in the case of long-range interactions. 

We now discuss the results from the simulations. Since there are strong arguments, 
discussed above, that z = 1, we work with systems with shape L x L x L with L = 6, 8, 10 
and 12. In most cases we perform 3,000 MCS for equilibration followed by 3,000 MCS with a 
measurement every 10 MCS. Close to the critical point the number of sweeps was generally 
larger for the larger sizes. The number of disorder realizations varied from 200 to 1000. As 
was the case for short-range interactions we carefully check for equilibration by computing 
the "Hamming distances" discussed in Section |V|. 



A. Determination of the Critical Point 



Since z = 1 it follows from the discussion after Eq. ( f4.5|) that we should look for the 
intersections of data for p(0)L against K for different lattice sizes. Our results are presented 
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in Fig. [13] for the case of the Ewald-sum form of the potential with e* = 1/2. Clearly 



the lines cross close to K — 0.240 and, more precisely, we estimate the critical coupling 
to be K c = 0.240 ± 0.003, quite close to the value for the short-range case. Since all four 
sizes intersect very close to the same point, Fig. [13| provides strong evidence that z — 1, in 
agreement with the scaling arguments. 

An equivalent analysis can be performed with the Green's function form for the potential, 
Eq. ( |2.44D , with e* 2 = 1/4. The two forms of the potential are different on short length scales, 
and the values of e* 2 are therefore not directly comparable. Although we don't have enough 
data for the largest size, 12 x 12 x 12, to perform a conclusive analysis, we can determine 
the critical coupling to be K c ~ 0.275, somewhat higher than for the Ewald form, with 
reasonable certainty. 



B. The Conductivity 



Following the approach used above for the short-range case, we plot, in Fig. |TJ, the 
resistivity, R^ in units of Rq, against frequency at the critical point, K c = 0.240. This 
data is for e* 2 = 1/2, with the Ewald method used to evaluate the Coulomb potential. The 
collapse of the data is excellent, without any correction involving T/u n , which was used for 
the short-range case. This implies that we can interchange the two limits ui — > 0, and T — > 
as expected. As for the short-range case, the data varies linearly with u n , implying a Drude 
form for the conductivity. Making a least square fit to the data with u n /uj c < 0.44 we find 
R*u/Rq= 1.82(2). 

We can now try to investigate the universality of R^ by studying the model with e* 2 = 1/4 
and the potential evaluated by the Green's function method. Fig. ^ shows the resistivity 



at the critical point K c = 0.275, along with the data already presented in Fig. |14]. We see 
two interesting things. Firstly, the actual form of the resistivity as a function of frequency is 
clearly different in the two cases. However, the extrapolation to zero frequency is the same 
within the uncertainties. For the Green's function potential the best fit is R^/Rq = 1.91(7), 
again for points with abscissa less than 0.44. The agreement between the zero-frequency 
limit of the two sets of data in Fig. [TB| provides strong support for the resistivity being 
universal at the transition. 

We have also studied the effect of the aspect ratio on the resistivity. This is important 
because the aspect ratio is related to quantities relevant for experiments, as we shall now 
show. An important concept in mesoscopic physics is the phase coherence length £ inc . This 
length is expected to diverge as T — > 0, and so, at criticality, should be proportional to 
the Bose glass correlation length, £, making the usual assumption that there is only one 
divergent length scale. To determine how £ varies as a function of T ~ L^ 1 note that from 
finite-size scaling, the finite-size relaxation time £ T at criticality should be proportional to L T 
and so characteristic lengths should scale as ^ z . Consequently £ inc ~ T -1 / 2 , which means 
that the the aspect ratio can be expressed as 

c = h. „ J (S*EL\ (7.11) 

i.e. it is proportional to a power of the ratio of the phase coherence length to the lattice 
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size. Experiments are generally carried out in the range L ^> £j nc so one should view 
the conductivity as arising from incoherent self- averaging of domains whose size is £i nc . In 
the opposite limit, L <C £ inc , we expect large variations from sample to sample ("bosonic 
universal conductance fluctuations"), and the average conductance will not necessarily be 
the same as that obtained in the regime L 3> £i nc . Thus, the conductivity at zero frequency 
but finite temperature is given by a scaling function, a(l/TL z ), where the argument is 
proportional to the the aspect ratio. The experimental situation corresponds to the limit of 
zero aspect ratio, whereas, the simulations are done for a finite value. 

We have therefore performed calculations for two other aspect ratios, 1/2 and 3/2 re- 
spectively. In both cases we used the Ewald form of the potential with e* 2 = 1/2. In order 
to obtain scaling plots for aspect ratios different from 1 it is necessary to include corrections 
of the form 1/L 2 . Including this correction term, our estimates for agree with those for 
aspect ratio unity, within the errors. Thus any dependence of Rjj on aspect ratio seems to 
be quite small. 

In conclusion, we estimate the universal conductivity at the critical point from all our 
data to be: 

R* D = (1.82 ± 0.20)R Q , a* = (0.55 ± 0.06)G Q . (7.12) 

No dependence on the aspect ratio, the microscopic form of the potential, the strength of 
the Coulomb interaction, or particle density was observed. 



C. The Wave- Vector Dependent Compressibility 

Because of the long-range interactions, the system is incompressible. As a result the 
wave-vector dependent compressibility should vary at criticality as n(k) ~ k from Eq. ([3.18 ) 
with z = 1. Fig. shows the data at the critical point K c = 0.240, for aspect ratio 1, and 
the Ewald form of the potential with e* 2 = 1/2. The solid lines shown are cubic splines 
fitted to the data points. The data for small k seems to be roughly linear, as expected if 
z = 1 but one would need substantially smaller wave vectors to draw a firm conclusion. 

Similar results results were obtained for aspect ratios 1/2 and 3/2. 



D. The Exponent v 

To determine the correlation length exponent we try to collapse the data in a scaling 
plot of p(0)L versus SL 1 ^, based on Eq. (|4.4j). Fig. [T7| shows the data for v = 0.90 and 



K c = 0.240, with aspect ratio 1, and the Ewald form of the potential with e* 2 = 1/2. By 
considering all our data we estimate 

v = 0.9 ±0.15. (7.13) 

The estimate of v is again consistent with the inequality of Chayes et al.il, v > 2/d, being 
satisfied as an equality. 
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E. The Correlation Functions 



Assuming z = 1 we expect the spatial correlation functions to have an asymptotic form of 
C x (r) ~ r~ Vx with y x = 1 + rj, see Eq. (|4.21 ). We have calculated the equal-time correlation 



functions for three different system sizes, L x L x L with L — 8, 10 and 12. The Ewald form 
of the potential was used with e* 2 = 1/2, and the calculation was performed at the critical 
point K c = 0.240. Despite there being quite a lot of noise in the data at large arguments, 
we can fit to the form in Eq. ( |6.11| ) with the result y x = 1.8 ± 0.4. Assuming that z = 1 
this then tells us that r] = 0.8 ± 0.4 in agreement with the previously derived inequality, 
Eq. (^TUp 7] < az, with a > 1. 

In principle, the time-dependent correlation functions can also be determined. However, 
unlike the equal-time correlation functions, these involve injecting an extra particle and 
then destroying it at a later time. Long-range interactions complicate the simulation of this 
correlation function and we have not attempted it. Hence we cannot independently confirm 
the value z — 1 of the dynamical exponent found in the scaling of p(0)L z , but believe it to 
be accurate. Similarly, we have not attempted to obtain the time decay of the correlation 
functions deep in the Bose glass phase, which could give us information on the Coulomb gap 
exponent, a in Eq. (|7.6| ). 



VIII. DISCUSSION 

We have investigated universal properties of the T = Bose glass to superfluid transition 
in two dimensions both for particles with short-range interactions and for particles with long- 
range (1/r) Coulomb interactions. We used a version of the path integral approach which 
corresponds to including only phase fluctuations of the condensate. 

For the case of long-range Coulomb interactions we find: v = 0.90 ± 0.15, rj = 0.8 ± 
0.4, z ~ 1.0, and a* = (0.55 ± 0.06) G Q where Gq 1 = Rq = h/(2e) 2 « 6.45 kfi. This model 
should be in the right universality class to describe the superconductor-insulator transition 
in disordered thin films. Experimental result c3 do not show much support for the 

universality of a*, the values of a* /Gq varying from about 0.6 to 2. However, it is possible 
that many of these experiments were not at sufficiently low temperature to probe the critical 
regime. We are not aware of any other theoretical calculations with which to compare our 
results. 

For the case of short-range interactions we find v = 0.9 ± 0.1, rj = —0.10 ± 0.15, and 
z = 2.0 ± 0.10. The universal conductivity is in this case a* = (0.14 ± 0.03)Gq. These 
results are in reasonable accord with those of RungoU, who studied a hard-core boson 
model on small lattices of size up to 4 x 4 by diagonalization, and obtained v = 1.4 ± 0.3 
and z = 1.95 ± 0.25. He also found a* /Gq = 0.16 ± 0.01, in agreement with our estimate, 
though the error bar may be rather optimistic, since different assumptions for the finite-size 
corrections led to significantly different values. 

Since this work was completed, two groups have reported results for the short-range 
model which differ from ours. Batrouni et alE3 have performed world-line quantum Monte 
Carlo simulations directly on the boson Hubbard model. They find er* = (0.45±0.07) Gq, but 
do not report values for the exponents. We do not have an explanation for this discrepancy, 
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though we note that their data for L 2 p(0) do not splay out in the insulating phase, K < K c , 
as do ours (see Fig. [^), which makes the location of the critical point harder. 

Results which are totally different from ours and also different from those of Batrouni et 
al.@ have been found by Makivic et al.H3 who performed world-line quantum Monte Carlo 
simulations on a hard-core Bose system finding er* = (1.2 ± 0.2) Gq, with z = 0.5 ± 0.05 
and v = 2.2 ± 0.2. They used a large lattice of size 64 x 64 x 48 in the zero winding number 
sector, and did finite-size scaling by considering sub-regions of sizes L su b x L su b x 48 with 
L su b between 4 and 32. Now the equilibration time varies as L 2mc , where, for this model, 
the Monte Carlo dynamical exponent, z mc , isill ~ 6. It is therefore surprising to us that 
such a large lattice can be equilibrated in the 3 x 10 5 sweeps that were used. It is also 
unclear if their data would not have been consistent with a larger dynamical exponent, 
had they scaled the the size of the sub-regions in the time direction like L* uh (and taken 
different values for z), rather than leaving this size fixed at L T = 48. Furthermore, the two 
temperatures used seem to be rather high since the critical coupling changed by a factor 
of two between them. Makivic et al.ill propose that the difference between their results 
and ours is that amplitude fluctuations, neglected in our model but included in the boson 
Hamiltonians, are relevant, so the two models are in different universality classes. While 
this idea certainly can not be ruled out, we don't yet feel that it has been conclusively 
demonstrated. First of all, the only evidence for it is the numerical results, about which we 
have some reservations discussed above. It would be more compelling if there were additional 
evidence, such as a calculation of the exponent for amplitude fluctuations, showing that it 
is indeed relevant in the renormalization group sense. Furthermore, this explanation does 
not explain the differences between the results of Makivic et al.E3 and those of Batrouni et 
and Rungeta, which also included amplitude fluctuations. 

For the future, it would be very interesting to study the field-tuned transitioning by 
Monte Carlo simulations, since this is expected to be in a different universality classS 
from the disorder-tuned transition discussed here. The problem is that one needs to find 
a representation of the problem in terms of a real classical (D+l) dimensional effective 
Hamiltonian which incorporates both the magnetic field and lack of microscopic particle-hole 
symmetry. Unfortunately, the phase representation, Eq. (|2.27| ), though it can be generalized 
to include a field and is still real for the particle-hole symmetric case, is complex in the 
absence of particle-hole symmetry, and the link representation used here becomes complex 
in the presence of a field, as can be seen from Eq. ( 2.35 ). We are therefore unaware of any 
representation of the problem suitable for Monte Carlo simulations. 
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FIGURES 




FIG. 1. Sketch of the general phase diagram for thin- film superconductors as a function of 
disorder, A, and temperature, T. The dashed line, T c q, indicates the mean-field onset temperature 
where Cooper pairs start to form. At T = an insulating phase appears (solid line) and a transition 
from a superconductor to an insulator takes place at a critical value of the disorder A c . 
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FIG. 2. Typical closed loop of integer-values currents on the links of the space-time lattice. 
A particle hops from x\ to X2 at time t±, diffuses until time T2 when it annihilates the hole it left 
behind. 
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FIG. 3. Schematic picture of a local move (lower left corner), and global move (right side of 
the picture). In the local moves changes of ±1 are attempted in the currents circulating around 
a plaquette. In the global the current is changed by ±1 along straight lines across the system, 
thereby changing the winding number, if the line is along a space direction, or the boson number if 
the line is in the time direction. The numbers indicate how the link variables are changed. Many 
local moves will change the current around an arbitrary loop as indicated in the upper left corner. 
A global move which destroys a boson is shown, and also another global move which increases the 
winding number along the x direction by one. 
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FIG. 4. The x and r Hamming distances for a system of size 8 x 8 x 16, with short-range 
interactions at the critical coupling K = 0.248. For each pair of curves, the upper one is for h a ^ 
and the lower one for h a . 
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FIG. 5. p(0)L 2 for the system sizes indicated, as a function of K for short-range interactions. 
From the intersection of the curves we estimate the critical coupling to be be K c = 0.248 ± 0.002 
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FIG. 6. The compressibility at zero wave vector, k(0), for different system sizes, as a function 
of K for short-range interactions. The critical point is at K c = 0.248. 
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FIG. 7. The compressibility of the model with short-range interactions as a function of wave 
vector, at the critical point, K c = 0.248, for different system sizes. The solid lines are spline fits to 
the data points, and k c = 2tt. 
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FIG. 8. The correlation function (imaginary time Green's function) for positive (imaginary) 
times, corresponding to a boson propagating forward in time, for a system of size 8 x 8 x 16 with 
short-range interactions, at K = 0.175 which is far into the Bose glass phase. The dashed line 
indicates a fit to the data of the form 0.170(2) (i-- 1 - 10 ® + (L T - t)- l10 ( 8 )). The data for r > L T /2 
is actually the value of C__(L T /2 — t) as discussed in the text. Thus the Green's function decays 
with a power of time (rather than exponentially) in the Bose glass phase, as predicted in Ref. [l0|. 
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FIG. 9. The resistivity in units of Rq = h/(2e) 2 , as a function of uj n /uj c — aL~ 2 (oj c /uj n ) 
for the model with short-range interactions, where to c = 2tt and a = 0.179. The calculation was 
done at the critical point, K = 0.248. The aspect ratio was in this case 1/4, and the system sizes 
shown were as indicated in the figure. The dashed line indicates a least square fit to the points 
with abscissa less than 0.26 of the following form 7.84(7) + 34.9(5)(u; ri /u; c — aL~ 2 (uj c /u n )). The 
correction involving the parameter a is proportional to T/u n as discussed in the text. On physical 
grounds we expect that a = in the thermodynamic limit, and indeed a fairly good fit, with almost 
the same value of the d.c. resistivity, is obtained with a = 0. 
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FIG. 10. Scaling plot of p(0)L 2 versus bL x l v , for short-range interactions, where S is the 
reduced coupling constant (K — K c )/K c and L is the linear system size. The parameters used in 
the plot are K c = 0.248, v = 0.90. 



47 



1.00 



"i r 



□ 10x10x25 
. 8x8x16 



"i i i r 



0.10 



0.01 



/ / 
/ / 

i 4 

/ i 
/ i 

■ m 

/ 



J I I I I I L 



10 



FIG. 11. The equal-time correlation function of the x-link variables as a function of spatial 
distance r(= x) for short-range interactions. The error bars increase with increasing r because the 
calculation involves the average of the exponential of a "string" of link variables, see Eq. ( |4.13 ), 
which can fluctuate hugely when the string is long. For this reason, the data for r > L/2 were deter- 
mined from C X (L — r). Two data sets are shown, for systems of size L = 8 and 10 with aspect ratio 
1/4, at K = 0.248, which is the critical point. The dashed line indicates fits to the data. For L = 8 
the fit is 0.18(1) (r- 2 - 02 W + {L — r)" 2 - 02 ^)). For L = 10 the fit is 0.18(l)(r- L94 ( 2 ) + (L - r)" L94 ( 2 )). 
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FIG. 12. The correlation function (Green's function) for positive times for short-range inter- 
actions for systems of size 6x6x9 and 8 x 8 x 16, i.e. with aspect ratio 1/4. The calculation was 
performed at the critical point, K c = 0.248. The dashed line indicates fits to the data. For L = 6 
the fit is 0.266(4)(r- L03 ( 1 ) +(L-r)- L03 ( 1 )). For L = 8 the fit is 0.256(4)(r-°- 94 ( 1 ) + (L-r)" - 94 ^)). 
The data for r > L T /2 is actually the value of C_(L r /2 — r), as discussed in the text. 
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FIG. 13. p{0)L for the system sizes indicated, as a function of K. The system sizes indi- 
cated correspond to an aspect ratio of 1. The critical point, determined from the intersections, is 
K c = 0.240 ± 0.003. The Ewald-sum form of the Coulomb potential was used with e* 2 = 1/2. 
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FIG. 14. The resistivity in units of Rq = h/(2e) 2 , as a function of io n /uj c , where uj c = 2tt. The 
calculation was done at the critical point, K = 0.240. The aspect ratio was in this case 1, and the 
system sizes shown were as indicated in the figure, e* 2 = 1/2 was used along with the Ewald form 
for the potential. The dashed line indicates a least square fit to the points with abscissa less than 
0.44 of the form 1.82(2) + 17.84(7)w„/u; c . 
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FIG. 15. The resistivity in units of Rq = /i/(2e) 2 , as a function of u! n /io c , where lo c = 2tt. 
For the upper curve the Green's-function form of the potential was used with e* 2 = 1/4, and an 
aspect ratio of 1. The calculation was done at the critical point for this potential, K = 0.275. 
The dashed line indicates a least square fit to the points with abscissa less than 0.44 of the form 
1.91(7) + 21.5(3)w n /w c . The lower curve is the results from the Ewald form of the potential from 
Fig. [14|. Although the results for the two forms of the potential differ at finite frequency, they 
appear to extrapolate to the same value in the d.c. limit, as expected since the d.c. resistivity is 
predicted to be universal. 
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FIG. 16. The compressibility at the critical point for the same model as in Fig. 14 



as a 



function of wave vector, for different system sizes, where k c = 2ir is the lattice cutoff. As expected 
in Coulomb systems, the compressibility appears to vanish as k — > 0. The solid lines are spline 
fits to the data points, and k c = 2tt. Also indicated by a dashed line is a tentative fit to the low 
frequency part with abscissa less than 0.19 of the form —0.020(2) + 0.8(2)k/k c . 
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FIG. 17. Scaling plot of p(0)L versus where 5 is the reduced coupling constant 

(K — K c )/K c and L is the linear system size. The model is the same as in Fig. [14]. The pa- 
rameters used in the plot are K c = 0.240, v = 0.90. 



54 



